Preliminary tests for dwarf galaxies orbits' plots with plotly
(given their ICRS coordinates)

from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import QTable
import astropy as a0
import numpy as np
import matplotlib.pyplot as plt
import math
import plotly
import plotly.graph_objs as go
import plotly.express as px

from visorbits import rotations as rot
from visorbits import ellipse as consec
plotly.offline.init_notebook_mode()

1. Initialize coordsiand transform to Galactocentric Cartesian system

---

c = {'Draco': {'ra': 15.*(17+20/60.+ 12.4/3600.)*u.deg, 'dec': (57+54/60.+55/3600.)*u.deg,\
                    'D':  75.8*u.kpc, 'v_r': -291.0*u.km/u.s,\
                    'mu_ra_cosd':   44.*10**(-3)*u.mas/u.yr, 'mu_dec': -188.*10**(-3)*u.mas/u.yr},\
       'Sculptor': {'ra': 15.*(1+  0/60.+  9.4/3600.)*u.deg, 'dec': (-33-42/60.-33/3600.)*u.deg,\
                    'D':  83.9*u.kpc, 'v_r':  111.4*u.km/u.s,\
                    'mu_ra_cosd':  100.*10**(-3)*u.mas/u.yr, 'mu_dec': -158.*10**(-3)*u.mas/u.yr},\
         'Fornax': {'ra': 15.*(2+ 39/60.+ 59.3/3600.)*u.deg, 'dec': (-34-26/60.-57/3600.)*u.deg, \
                    'D': 147.2*u.kpc, 'v_r':   55.3*u.km/u.s,\
                    'mu_ra_cosd':  381.*10**(-3)*u.mas/u.yr, 'mu_dec': -359.*10**(-3)*u.mas/u.yr},\
          'LeoII': {'ra': 15.*(11+13./60.+28.8/3600.)*u.deg, 'dec': (22+ 9/60.+ 6/3600.)*u.deg, \
                    'D': 233.0*u.kpc, 'v_r':   78.0*u.km/u.s,\
                    'mu_ra_cosd': -109.*10**(-3)*u.mas/u.yr, 'mu_dec': -150.*10**(-3)*u.mas/u.yr},\
           'LeoI': {'ra': 15.*(10+ 8./60.+28.1/3600.)*u.deg, 'dec': (12+18/60.+23/3600.)*u.deg,\
                    'D': 258.2*u.kpc, 'v_r':  282.5*u.km/u.s,\
                    'mu_ra_cosd':  -50.*10**(-3)*u.mas/u.yr, 'mu_dec': -120.*10**(-3)*u.mas/u.yr}}

set_of_icrs = {}

for dSph in c:
    set_of_icrs[dSph] = SkyCoord(ra=[c[dSph]['ra']], dec=[c[dSph]['dec']], distance=[c[dSph]['D']],\
                                 pm_ra_cosdec=[c[dSph]['mu_ra_cosd']], pm_dec = [c[dSph]['mu_dec']],\
                                 radial_velocity=[c[dSph]['v_r']], frame='icrs')
    
    sGC0 = set_of_icrs[dSph].transform_to(a0.coordinates.Galactocentric)

    if dSph=='Draco':
        sGC = sGC0.to_table()
        sGC = a0.table.hstack([sGC, QTable({'Name': [dSph]}, names=['Name'])])
    else:
        sGC = a0.table.vstack([sGC, a0.table.hstack([sGC0.to_table(), QTable({'Name': [dSph]}, names=['Name'])])], join_type='inner')

data = go.Scatter3d(x=sGC['x'], y=sGC['y'], z=sGC['z'], text = sGC['Name'], mode='markers', marker=dict(size=4), name='Dwarves')

fig = go.Figure(data)

for i, n in enumerate(sGC['Name']):
    fig.add_trace(
            go.Scatter3d(
                x=[sGC['x'][i].value, sGC['x'][i].value+sGC['v_x'][i].value],
                y=[sGC['y'][i].value, sGC['y'][i].value+sGC['v_y'][i].value],
                z=[sGC['z'][i].value, sGC['z'][i].value+sGC['v_z'][i].value],
                mode='lines', legendgroup=n,  hovertext=n,
                name="{}".format(n)  )  )

fig.update_layout(title='Figure 1.', 
                  autosize=True,
                  scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
                  scene = dict(
                    xaxis_title='X', yaxis_title='Y', zaxis_title='Z',
                    xaxis = dict(nticks=15,range=[-770,770],),
                    yaxis = dict(nticks=15,range=[-770,770],),
                    zaxis = dict(nticks=15,range=[-770,770],),),
                  width=550, height=550, template='none', scene_aspectmode='cube')
fig.show()

2. Prepare plotly local functions

---

Function	Arguments	Output	Description
psurface	point p0 (3,) point p1 (3,)	Vector in the place, tau Normal vector X, Y of points (3,) ZV (meshgrid of plane z-points for xiven X, Y)	Define our plane which contains the zero point
rotatedset	initial (3, len(p)) set point p1 (3,) tau_ort (3,) vector which lies in the plane (perpendicular to p1 as example) phi4 -- phase angle for final rotation	initial ell1, ell2, ell3, ell4 -- steps of rotations r1, r2, r3, r4 -- rotational matrices	Step-bystep rotation of initil dataset.
initfig_with_surface	-- (len(p),) array of x y -- (len(p),) array of y zv -- meshgrid array of z	fig object of plotly	Initilizes fig
plotvectors	fig (plotly object; I'm not familiar with it) vectors (array of pairs of (3,) points for beginning and end of each vector) vc_legends (list of strings for legends) vc_colors (list of colors)	--	Updates fig
plottraces	fig, traces, sys_legends, colors fig (plotly object; I'm not familiar with it) traces (list of all initial and rotated lines, (len(p), 3)) sys_legends (list of strings for legends) colors (list of colors)	--	Updates fig

def psurface(p0, p1):
    # p0 and p1 together with (0,0,0) define the plane.
    # Find the normal vector:
    norm = rot.vect_m(p0, p1)
    norm /= np.sqrt(np.sum(norm**2))
    # vtau_test is the direction of motion
    # p1 marks the location of our point (which lies on the future trace)
    vtau_test = rot.vect_m(norm, p1)
    #tau_ort  = rot.vect_m(norm_ort, r_ort)  # vtau_test
    
    tau = vtau_test/np.sqrt(np.sum(vtau_test**2))
    # Definition of the final XY plane:
    x = np.array([0, p0[0], p1[0]])
    y = np.array([0, p0[1], p1[1]])
    
    z0      = lambda xl, yl: (xl*norm[0] +yl*norm[1])*(-1)/norm[2]
    xv, yv  = np.meshgrid(x,y)
    zv      = z0(xv, yv)
    return norm, tau, x, y, zv

def rotdataset(initial, p1, tau_ort, phi4):
    # Create roational matrices 
    r1, r2, r3, r4 = rot.reverse_tr(p1[0], p1[1], p1[2], tau_ort, phi4=phi4)   
    # Application of rotational matrices, step-by-step:
    ell1 = np.dot(r1, initial.T).T
    ell2 = np.dot(r1, np.dot(r2, initial.T)).T
    ell3 = np.dot(r1, np.dot(r2, np.dot(r3, initial.T))).T
    ell4 = np.dot(np.dot(r1, np.dot(r2, np.dot(r3,r4))), initial.T).T
    return initial, ell1, ell2, ell3, ell4, r1, r2, r3, r4

def initifig_with_surface(x, y, zv):
    # Figure: vectors and 3-axes of coordinates.
    fig = go.Figure(data=[go.Surface(x=x, y=y, z=zv,\
                                     surfacecolor=np.full(zv.shape,'green'),opacity=0.2, showscale = False, name='Destination')])
    return fig

def plotvectors(fig, vectors, vc_legends, vc_colors):
    for i, v in enumerate(vectors):
        fig.add_trace(go.Scatter3d(x=[v[0][0], (v[1][0]-v[0][0])*3],\
                                   y=[v[0][1], (v[1][1]-v[0][1])*3],\
                                   z=[v[0][2], (v[1][2]-v[0][2])*3], mode='lines',\
                                   line = dict(color=vc_colors[i]) ,line_width=3,showlegend=True,name=vc_legends[i],))
def plottraces(fig, traces, sys_legends, colors):
    for i, d in enumerate(traces):
        fig.add_trace(go.Scatter3d(x=d[:,0], y=d[:,1], z=d[:,2],\
                                   line = dict(color=colors[i]),\
                                   marker_size=1,showlegend=True,name=sys_legends[i],))

3. Tests.

---

3.1 Three-axes:

# Initialize two points, which with (0,0,0) determine the plane of final XY
p0x =  0.201; p0y  = -0.8;        
p1x = -0.1;   p1y = 0.85

p0  = np.array([p0x, p0y, np.sqrt(1-p0x**2-p0y**2)])
p1  = np.array([p1x, p1y, np.sqrt(1-p1x**2-p1y**2)])

initial = np.array([[0,1,0,0,0,0],\
                    [0,0,0,1,0,0],\
                    [0,0,0,0,0,1]]).T

norm, tau, x, y, zv =  psurface(p0, p1)
initial, ell1, ell2, ell3, ell4, r1, r2, r3, r4 = rotdataset(initial, p1, tau, phi4=0)

fig0 = initifig_with_surface(x, y, zv)

# Prepare settings for vectors' plot:
vectors    = np.array([[np.zeros(3,), norm], [np.zeros(3,), p1], [p1, p1+tau]])
vc_legends = ['Normal', 'Rad-vector of the point which belongs to the final ellipse (x3)', 'Tangential velocity direction']
vc_colors  = ['pink', 'orange', 'cyan']

# Prepare settings for lines' plot:
traces = list([initial, ell1, ell2, ell3, ell4])
sys_legends = ['Initial', 'Aligned with XY-projection of point, which belongs to final',\
               'Final X is he point radius-vector', 'XY is in the correct plane', 'Correct phase']
colors      = ['black', 'red', 'orange', 'lime', 'cyan']
  

plotvectors(fig0, vectors, vc_legends, vc_colors)
plottraces(fig0, traces, sys_legends, colors)

fig0.update_layout(title='Figure 2. Test rotation of 3-axes.', autosize=True,
                  scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
                  scene = dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z',\
                               xaxis = dict(nticks=5), yaxis = dict(nticks=5), zaxis = dict(nticks=5),),\
                  width=1000, height=700, template='none', scene_aspectmode='data')
fig0.show()

3.2 Simple ellipse:

# Initialize two points, which with (0,0,0) determine the plane of final XY
p0x =  0.201; p0y  = -0.8;        
p1x = -0.1;   p1y = 0.85

p0  = np.array([p0x, p0y, np.sqrt(1-p0x**2-p0y**2)])
p1  = np.array([p1x, p1y, np.sqrt(1-p1x**2-p1y**2)])

norm, tau, x, y, zv =  psurface(p0, p1)

# Definition of the orbit in its own space:
# pv0 - is our curent point with phi4 phase angle from pericenter
ang = np.linspace(-np.pi, np.pi, 50); 
e   = 0.6
r = np.sqrt(np.sum(p1**2))/(1-e)
xcircle = r*(1-e**2)/(1+e*np.cos(ang))*np.cos(ang)
ycircle = r*(1-e**2)/(1+e*np.cos(ang))*np.sin(ang)
zcircle = np.zeros(len(ang))

points1 = np.hstack([np.reshape(xcircle, (len(xcircle), 1)), np.reshape(ycircle, (len(xcircle), 1))])
initial = np.hstack([points1, np.reshape(zcircle, (len(xcircle), 1))])

initial, ell1, ell2, ell3, ell4, r1, r2, r3, r4 = rotdataset(initial, p1, tau, phi4=0)

fig1 = initifig_with_surface(x, y, zv)

# Prepare settings for vectors' plot:
vectors    = np.array([[np.zeros(3,), norm], [np.zeros(3,), p1], [p1, p1+tau]])
vc_legends = ['Normal', 'Rad-vector of the point which belongs to the final ellipse (x3)', 'Tangential velocity direction']
vc_colors  = ['pink', 'orange', 'cyan']

# Prepare settings for lines' plot:
traces = list([initial, ell1, ell2, ell3, ell4])
sys_legends = ['Initial', 'Aligned with XY-projection of point, which belongs to final',\
               'Final X is he point radius-vector', 'XY is in the correct plane', 'Correct phase']
colors      = ['black', 'red', 'orange', 'lime', 'cyan']

plotvectors(fig1, vectors, vc_legends, vc_colors)
plottraces(fig1, traces, sys_legends, colors)

fig1.update_layout(title='Figure 3. Test with ellipses.', autosize=True,
                  scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
                  scene = dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z',\
                               xaxis = dict(nticks=5), yaxis = dict(nticks=5), zaxis = dict(nticks=5),),\
                  width=1000, height=700, template='none', scene_aspectmode='data')
fig1.show()

3.3 All rotations of our data:

s = sGC

s = sGC
ellipses = {}
for key in list(sGC['Name']):
    ellipses[key] = {'x': [], 'y': [], 'z': [], 'set': []}

for i00, n in enumerate(s['Name']):
    
    r_dwarf = np.array([s['x'][i00].value, s['y'][i00].value, s['z'][i00].value])
    v_dwarf = np.array([s['v_x'][i00].value, s['v_y'][i00].value, s['v_z'][i00].value])
    
    r0 = np.sum(r_dwarf**2)**0.5;
    v0 = np.sum(v_dwarf**2)**0.5
    
    r_ort    = r_dwarf/r0
    #norm_ort = rot.vect_m(r_ort, v_dwarf)/np.sum(rot.vect_m(r_ort, v_dwarf)**2)**0.5
    norm_ort, tau_ort_wrong, x, y, zv =  psurface(r_dwarf, v_dwarf)
    tau_ort  = rot.vect_m(norm_ort, r_ort)  # vtau_test
    
    vr = v0*np.sum(r_ort*v_dwarf)/(np.sum(r_ort**2)**0.5*np.sum(v_dwarf**2)**0.5) * r_ort
    vt = tau_ort * np.sqrt(v0**2 - np.sum(vr**2))
    vt0 = np.sqrt(np.sum(vt**2))
    
    e = consec.ecc(r0*u.kpc, v0*u.km/u.s, vt0*u.km/u.s)
    a = consec.a_ax(r0*u.kpc, v0*u.km/u.s)
    
    psi_current = consec.angl_ecc(r0, e, a)
    while abs(psi_current)>2*np.pi:
        psi_current +=2*np.pi*(-1*psi_current)/(abs(psi_current))
    
    phi4 = consec.angl_c(psi_current, e)
    while abs(phi4)>2*np.pi:
        phi4 +=2*np.pi*(-1*phi4)/(abs(phi4))

    # Definition of the orbit in its own space:
    # r_dwarf - is our curent point with phi4 phase angle from pericenter
    ang     = np.linspace(0, np.pi*2, 50);
    xcircle = a*(1-e**2)/(1+e*np.cos(ang))*np.cos(ang)
    ycircle = a*(1-e**2)/(1+e*np.cos(ang))*np.sin(ang)
    zcircle = np.zeros(len(ang))
    points1 = np.hstack([np.reshape(xcircle, (len(xcircle), 1)), np.reshape(ycircle, (len(xcircle), 1))])
    initial = np.hstack([points1, np.reshape(zcircle, (len(xcircle), 1))])

    p1 = np.array([xcircle[0], ycircle[0], zcircle[0]])
    p2 = np.array([a*(1-e**2)/(1+e*np.cos(phi4))*np.cos(phi4),\
                  a*(1-e**2)/(1+e*np.cos(phi4))*np.sin(phi4), 0.])
    
    # Create roational matrices 
    initial, ell1, ell2, ell3, ell4, r1, r2, r3, r4 = rotdataset(initial, r_dwarf, vt, phi4=phi4)
    
    ellipses[n]['x'] = ell4[:,0];   ellipses[n]['y'] = ell4[:,1];     ellipses[n]['z'] = ell4[:,2]; ellipses[n]['set'] = ell4

    p1 = np.dot(np.dot(r1, np.dot(r2, np.dot(r3,r4))), p1)
    p2 = np.dot(np.dot(r1, np.dot(r2, np.dot(r3,r4))), p2)

    fig2 = go.Figure()
    
    # Prepare settings for vectors' plot:
    vectors    = np.array([[np.zeros(3,), norm_ort*a], [np.zeros(3,), r_dwarf], [r_dwarf, r_dwarf+v_dwarf],\
                           [np.zeros(3,), p1], [np.zeros(3,), p2]])
    vc_legends = ['Normal', 'Rad-vector of the point which belongs to the final ellipse (x3)',\
                  f'Velocity direction, phase angle is {phi4*180./np.pi}', 'Pericenter', 'The Dwarf']
    vc_colors  = ['pink', 'orange', 'cyan', 'black', 'red']
    
    # Prepare settings for lines' plot:
    traces = list([initial, ell1, ell2, ell3, ell4])
    ell_legends = ['Initial', 'Aligned with XY-projection of point, which belongs to final',\
                   'Final X is he point radius-vector', 'XY is in the correct plane', 'Correct phase: X is in the pericenter']
    colors      = ['black', 'red', 'orange', 'lime', 'cyan']

    plotvectors(fig2, vectors, vc_legends, vc_colors)
    plottraces(fig2, traces, ell_legends, colors)
    
    fig2.update_layout(title=f'Figure {4+i00}. {n}', autosize=True,
                      scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
                      scene = dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z',\
                                   xaxis = dict(nticks=5), yaxis = dict(nticks=5), zaxis = dict(nticks=5),),\
                      width=1000, height=700, template='none', scene_aspectmode='data')
    
    fig2.show()

3.4 Plot everything:

fig3 = go.Figure()
colors_dwarf = ['red', 'orange', 'green', 'blue', 'purple']
for i, n in enumerate(list(sGC['Name'])):
    d = ellipses[n]['set']
    plottraces(fig3, list([d]), [n], [colors_dwarf[i]])
    fig3.add_trace(
            go.Scatter3d(
                x=[sGC['x'][i].value, sGC['x'][i].value+sGC['v_x'][i].value],
                y=[sGC['y'][i].value, sGC['y'][i].value+sGC['v_y'][i].value],
                z=[sGC['z'][i].value, sGC['z'][i].value+sGC['v_z'][i].value],
                mode='lines', legendgroup=n,  hovertext=n, line=dict(color=colors_dwarf[i]),
                name="{} position".format(n)  )  )

fig3.update_layout(title='Figure..', autosize=True,
                  scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
                  scene = dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z',\
                               xaxis = dict(nticks=5), yaxis = dict(nticks=5), zaxis = dict(nticks=5),),\
                  width=1000, height=700, template='none', scene_aspectmode='data')
fig3.show()